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Abstract: Subdivision surfaces are proven to be a powerful tool in geometric modeling and 
computer graphics, due to the great flexibility they offer in capturing irregular topologies. This 
paper discusses the robust and efficient implementation of an isogeometric discretization approach 
to partial differential equations on surfaces using subdivision methodology. Elliptic equations with 
the Laplace-Beltrami and the surface bi-Laplacian operator as well as the associated eigenvalue 
problems are considered. Thereby, efficiency relies on the proper choice of a numerical quadrature 
scheme which preserves the expected higher order consistency. A particular emphasis is on the 
robustness of the approach in the vicinity of extraordinary vertices. In this paper, the focus is 
on Loop’s subdivision scheme on triangular meshes. Based on a series of numerical experiments, 
different quadrature schemes are compared and a mid-edge quadrature, which is easy-to-implement 
via lookup tables, turns out to be a preferable choice due to its robustness and efficiency. 

Keywords: subdivision methods, isogeometric analysis, PDEs on surfaces, numerical quadrature. 

1 Introduction 

During the last years, isogeometric analysis (IgA) [U |2] emerged as a means of unification of the 
previously disjoint technologies of geometric design and numerical simulation. The former technology 
is classically based on parametric surface representations, while the latter discipline relies on finite 
element approximations. The advantages of the isogeometric paradigm include the elimination of 
both the geometry approximation error and the computation cost of changing the representation 
between design and analysis. A particular challenge is the problem of efficient numerical integration. 
Indeed, the higher degree and multi-element support of the basis functions seriously affect the cost of 
robust numerical integration. Consequently, numerical quadrature in IgA is an active area of research 

mmumm- 

Subdivision schemes are widespread in geometry processing and computer graphics. Eor a com¬ 
prehensive introduction to subdivision methods in general we refer the reader to [8] and [9]. With 
respect to the use of subdivision methods in animation see |10[ m- Nowadays, subdivision finite 
elements are also extensively used in engineering [laiiaiiiiiisiiiniiiTi. This way, the subdivision 
methodology changed the modeling paradigms in computer-aided design (CAD) systems over the 
last years tremendously, e.g. the Freestyle module of PTC Creo®, the Power Surfacing add-on for 
SolidWorks® or NX RealizeShape of Siemens PLM®. For the integration of subdivision methods 
in CAD systems, we refer to [lailg] and the references therein. With this development, powerful 
computer-aided engineering (CAE) codes have to be implemented for subdivision surfaces. Today 
NURBS and subdivision surfaces co-exist in hybrid systems. 

Among the most popular subdivision schemes are the Catmull-Clark [20) and Doo-Sabin |21j 
schemes on quadrilateral meshes, and Loop’s scheme on triangular meshes [22]. The Loop subdivision 
basis functions have been extensively studied, most interestingly regarding their smoothness |23[ I24| . 
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(local) linear independence [25l[26], approximation power [27], and the robust evaluation around so 
called extraordinary vertices |28| . 

Computing PDEs on surfaces is an indispensable tool either to approximate physical processes on 
the surfaces [29] or to process textures on the surface or the surface itself [30l |3ll [32] . Recently, surface 
PDEs have been considered in the field of IgA, and both numerical and theoretical advancements 
are reported [MllMlES]. 

The use of Loop subdivision surfaces as a finite element discretization technique is discussed in 
m and isogeometric discretizations based on the Catmull-Clark scheme in [MIET], whereas [38] 
investigates the use of the Doo-Sabin scheme in the finite element context. An isogeometric finite 
element analysis based on Catmull-Clark subdivision solids can be found in [39]. 



Figure 1: As an application we show here the first 24 eigenfunctions of the Laplace-Beltrami operator 
computed on a complex surface using the isogeometric subdivision approach. See Figure for the 
coarse control mesh. 

In this paper, we focus on the Loop subdivision scheme [ID] on triangular meshes. Loop subdivi¬ 
sion is known to describe limit surfaces A4 except at finitely many points (called extraordinary 
vertices) where they are only n [231124] . This allows us to use them in a conforming finite ele¬ 
ment approach not only for second-order but also for fourth-order elliptic PDEs. As model problem 
we consider the Laplace-Beltrami equation 

-AMU = fonM, ( 1 ) 

the surface bi-Laplacian equation 

{-AMfu = fonM ( 2 ) 

as well as the eigenvalue problem 

—Aj^u = \u on M. . (3) 

Let us remark that on closed, smooth surfaces the kernel of the Laplace-Beltrami and the surface 
bi-Laplacian operators are the constant functions. 

The crucial step in the implementation of a subdivision finite element method is the choice of 
numerical quadrature. The computational cost, consistency, robustness and the observed order of 
convergence are important parameters to evaluate the appropriateness of each numerical integra¬ 
tion technique. We aim at a detailed experimental study of the convergence behavior for different 
quadrature rules. We examine Gaussian quadrature, barycenter quadrature, mid-edge quadrature 
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and discuss adaptive strategies around extraordinary vertices. Furthermore, we provide a look-up 
table for the mid-edge quadrature that facilitates the implementation and results in a fairly robust 
simulation tool and a very efficient assembly of finite element matrices. This allows for a straightfor¬ 
ward extension of existing subdivision modeling codes to simulations with subdivision surfaces for 
applications, see Figure [T} 

Organisation of the paper. In Section we fix some notation and define subdivision functions 
and subdivision surfaces. The isogeometric subdivision method to solve elliptic partial differential 
equations on surfaces is described in Section where we derive the required finite element mass and 
stiffness matrices. Section investigates the different numerical quadrature rules and comments on 
the efficient implementation. In Section we report on the observed convergence behavior for a set 
of selected test cases. Finally, we draw conclusions in Section 

2 Subdivision Functions and Surfaces 

To discretize geometric differential operators and to solve geometric partial differential equations on 
closed subdivision surfaces we consider isogeometric function spaces defined via Loop’s subdivision. 

The domain of a Loop subdivision surface is a topological manifold obtained by gluing together 
copies of a standard triangle. More precisely, we consider the unit triangle A C with vertices 
(0,0), (1,0) and (0,1), which is considered as a closed subset of and a finite cell index setic C Z. 
The Cartesian product A x Xc defines the pre-manifold which consists of the cells A x {i}, i £ Xc. 
The manifold is constructed by identifying points on the boundaries of the cells, as described below. 

In addition to the cell index set, we consider an edge index setle dcX Xc, which is assumed to 
be symmetric 

(b j ) ^ ^ (jA) ^ 

and irreflexive 

\li ^TL\ {i, i) 0 Xg. 

Each edge is described by a symmetric pair of edge indices {{i,j), £ Zg x Xg. Moreover, each 

cell contributes to exactly three edges, 

Vi E Xc : \{j : {i,j) E Xg}] = 3. 



Figure 2: The spherical subdivision limit surface Spherical-5-12 with successively refined control 
meshes. The control meshes are plotted on the limit surface, i.e., they visualize the curved IgA- 
elements on the surface. Furthermore, the control meshes possess two vertices of valence 12 and 24 
vertices of valence 5 on every refinement level. Since two EVs share an edge in the coarsest control 
mesh, this mesh cannot be used for simulation. 

Each edge index {i,j) is associated with a displacement 6j, which is one of the 18 isometries 
(3 • 3 • 2) that map the unit triangle to one of its three neighbor triangles obtained by reflection 
and index permutation. These displacements in particular identify the common edges of neighboring 
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triangles and its orientation with respect to the two triangles. They need to satisfy the following two 
conditions. Firstly, the identification of the common edge and its orientation has to be consistent for 
each edge, thus 

Secondly, each edges of a triangle is identified with exactly one edge of another triangle, 
y{i,j) E le : y{i,k) E le : <5*(A) n A = Sl{A) nA^ j = k. 

Two points of the pre-manifold are identified, denoted by {^,i) ~ if (hj) is an edge index 

and the displacement (5®- transforms r] into This implies that ^ and rj are located on the boundary 
of the standard triangle. We denote with ~ the reflexive and transitive closure of this relation. This 
closure leads to the obvious identification of common vertices. The topological manifold (A x Tc)/~ 
is the domain manifold of the Loop subdivision surface. 

There, the indexed standard triangles define an initial triangulation (of level 0) of the domain 
manifold. Now, Loop’s subdivision scheme iteratively creates triangulations of level i > 0 containing 
4^|Xc| triangles. They are obtained by creating new vertices at the edge midpoints and splitting each 
triangle of level ^ — 1 into four smaller ones. The valence of the vertices in these triangulations are 
always 6, except for those vertices of the coarsest level that possess a different valence. The latter 
ones are called extraordinary vertices (EVs). For an example of a Loop subdivision surface with 
extraordinary vertices and different resolutions of the control mesh see Figure 

Next, we consider the spaces of piecewise linear functions on the triangulation of level £. 
Each function is uniquely described by its nodal values, i.e., by its values at the vertices of the 
triangulation. The Loop subdivision operators —)• are linear operators that transform 

a piecewise linear function of level £ — 1 into a function of level £. The nodal values of the 
latter (refined) function are computed as weighted averages of nodal values of the coarser function 
at neighboring vertices. Values at newly created vertices (at edge midpoints) depend on the four 
nodal values of the coarser function on the two triangles intersecting at the edge, while values at 
existing vertices are computed from the values at the (valence+l)-many vertices in the one-ring 
neighborhood. The weights are determined by simple formulas m- 

The space of Loop subdivision splines [40] consists of the limit functions generated by the subdi¬ 
vision operators, 

{ lim : / E L°}. 

L—¥ oo 

These limit functions are known to be bivariate quartic polynomials on all triangles of the triangu¬ 
lation of level i that do not possess any EVs. They are C^-smooth everywhere, except at the EVs. 
If one uses a special parameterization around EVs, which is given by the charateristic map [231 
then the limit functions are C^-smooth at EVs. 

Let ly be the vertex index set of the coarsest triangulation. Each vertex i ^Xy has an associated 
piecewise linear hat function Aj E L^ which takes the nodal value 1 at the associated vertex and 0 
else. The limit functions 

= lim R^R^-^ ■ ■ ■ R^Ki, 

L —^ oo 

which will be denoted as the Loop basis functions, span the space of Loop subdivision splines. The 
support of the Loop basis function is the two-ring neighborhood of the vertex i in the coarsest 
triangulation. If the support does not contain EVs in its interior, then it consists of 24 triangles 
of the coarsest level and the Loop basis function is the C^-smooth quartic box spline on the type- 
1 triangulation. Otherwise, the Loop basis function is still a piecewise polynomial function but 
possesses an infinite number of polynomial segments. These are associated with triangles obtained 
by successive refinement in the vicinity of the EV and not containing the EV itself [8]. 
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Now we have everything at hand to introduce a Loop subdivision surface 
X: Ax ^ m3 . k) = Y, Ci^C, k ), 

ieit, 

which is obtained by assigning control points Ci G M^ to the Loop basis functions with = (^i, .^ 2 ) £ 
A. This surface is smooth everywhere, except at EVs (in the same sense as described before) 
|23l I24j . Hence, it can be shown that the surface possesses a well-defined tangent plane at EVs, 
provided that the control points are not in a singular configuration. We denote the surface (as an 
embedded manifold) hy Ai = A4[V]. In the context of isogeometric analysis, the parameterization 
X of the domain is referred to as the geometry mapping. The basis functions (fi of the associated 
isogeometric function space are the push-forwards of the subdivision basis functions 

(Pi{x) =<PiO X~^{x), 


where x = X{^,k), i.e., the basis functions are defined by the discretization X. We define the 
discretization space as 

14 = G spanjgjjv 9 i}|. 

Here, h indicates the grid size of the control mesh, i.e., h = max(j \Ci — Cj\. For each cell, we 
denote by J{^,k) = (V^ X^ 2 ) the Jacobian of X and by G = J^J its first fundamental form, 
where X^i and X ^2 denote the tangent vectors with X^i = ^V(^, A;). Furthermore, for a function 
tt : AJ —M 


Vmu{x) = (JG-^Vg(no V)) {(,k) 
is the (embedded) tangential gradient or surface gradient and 

A_v(m(x) = divg ("Vdet GG~^'V^{u o X) \ {^,k) 

V det G ^ ' 

the Laplace-Beltrami operator, where x = X(^, k). 

3 Isogeometric Subdivision Method to solve PDEs 

Similar to the classical finite element method, the isogeometric subdivision approach transforms 
the strong formulation of a partial differential equation (e.g. 0-0) into a corresponding weak 
formulation on a suitable subspace of some Sobolev space and approximates the solution of the 
weak problem in the finite-dimensional sub-spaces C V^. We consider 


V^ = H\M)n{ uda = 0} 

Jm 

for problem ([^ and 

V^ = H^{M)n{[ uda = 0} 

Jm 

for problem Q. Following the general isogeometric paradigm 

v^ = Vhn{[ u;,da = 0} 

Jm 

as the discrete ansatz space on the subdivision surface At, e.g., see Figure [2j 


we use 
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Let a : x —>• M denote a symmetric, coercive and bounded bilinear form on and 

if :V^ a linear form for a given function / G Lp‘{Ai). Then there exists a unique solution 

u G of the variational problem 


a{u,v) = if{v), V u G 

where if{v) = fvda ( @21USJ HU |35] ). For our model problems, we use 


(4) 


and 


a{u,v) = / V_A 4 U • V_a 4 U da 

Jm 


a{u,v) = / Aji 4 uAj^vda 

Jm 


for the Laplace-Beltrami 0 and the Bi-Laplacian problem Q , respectively. The associated Galerkin 
approximation asks for the unique solution Uh G of the discrete variational problem 


a{uh,Vh) = ifivh), V u/j G V° 


(5) 


Due to the known Lf^-regularity of the Loop subdivision splines this Galerkin approximation is 

conforming, i.e., C V^. Using the basis expansion Uh = '^h with coefficients Ui £ R 

i&Xv 


one obtains the linear system 


SU = B, 


( 6 ) 


where U = {ui)i^x^ G denotes the coefficient vector, Sij = a{(pi,(pj) G the stiffness 

matrix, and Bj = Jj^fipjda G the right-hand side. The stiffness matrix for the Laplace- 

Beltrami problem Q is given by 

S*y= [ VMVi-^MV’jda=y [ V^^iG-^V^^.VdetGdC (7) 

and for the surface bi-Laplacian problem 0 we obtain 

Sg' = ^ (Vdet div^ (^Vdet GQ-^V 

The variational formulation of the eigenvalue problem 0 consists in finding a solution (u. A) G xM 
such that 


a{u, v) = X ■ m{u, v), V v G V^, 

where m{u,v) = uvda denotes the L^-product on A4. Furthermore, we denote by M the mass 
matrix with 

Mij = m{ipi,ipj) = y / ^i^jVdetGd,^. (9) 

feeic 

We obtain the discrete eigenvalue problem SU = A/^MU which can be solved by inverse vector 
iteration with projection @6] . 
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4 Numerical Quadrature 

In this section, we recap the evaluation-based assembly of the previously defined discrete variational 
formulations. We review standard Gaussian and barycentric quadrature assembly, discuss the special 
treatment at extraordinary vertices and finally introduce an edge-centered assembly strategy which 
is based on simple hxed, geometry-independent lookup tables and already very good consistency in 
our experiments (see Section]^. As discussed above, the isogeometric subdivision approach is based 
on higher order spline discretizations. The integrands of the corresponding IgA-matrices (Q ^ Q) 
are in general nonlinear and even on planar facets of the subdivision surface (with constant metric 
G) high-degree polynomials. Indeed, apart from EVs the limit functions d>i are quartic polynomials 
and thus for the mass matrix the integrand is a polynomial of degree 8 and for the stiffness matrix 
0 of degree 6 and for ([^ of degree 4. 

In the variational formulation, the integration is always pulled-back to the reference triangle, i.e., 

/ /da= V / /oAv/detG(e,A:)d^ 
fceic 


Thus numerical quadrature is applied to these pulled-back integrands. For a general function g using 
the evaluation-based quadrature 

K 

9=1 

with quadrature points on the reference triangle A and weights Wq, we replace the stiffness 
matrices, the mass matrix, and the right-hand side by the following quadrature based counterparts 



K 

Sg= ^ (10) 

feeic 9=1 


and 



div^ (^Vdet div^ (^Vdet ^ ) (C'', k) (11) 

fceXc 9 =i vdetG/ 


as well as 


K 

M.,= EE “’»(•*■ i-<PjVdetG) (^^A:) (12) 

keic 9=1 


B.= EE-.((/ O A) • Vdet g) (e^A:). (13) 

fceXc 9=1 

Now, instead of solving Q, we solve the system S • U = B. The associated, modified variational 
problem reads as follows: Find Uh £ such that 

a{uh,Vh) = ^f{vh), V Ufe G 14°, (14) 

where a{ipi, ipj) = Sjj and = Bj. If d(.,.) is uniformly 14-elliptic (S positive definite) existence 

and uniqueness of the discrete solution Uh is ensured (cf. [l2], for instance). 
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Gaussian Quadrature. 

For second-order elliptic problems, standard error estimates for finite element schemes with Gaussian 
quadrature |44l Chapter 4.1] imply that the expected order of convergence, in case of exact integra¬ 
tion, is preserved if the quadrature scheme is exact for polynomials of degree p = 6 . Transferring 
this to fourth-order problems we request exactness for polynomials of degree p = 4. This suggests to 
choose a Gaussian quadrature rule GA(p) which guarantees this required exactness. For the Laplace- 
Beltrami equation Q iF = 12 and for the surface bi-Laplacian equation Q iF = 6 quadrature points 
have to be taking into account on the reference triangle A (for symmetric Gaussian quadrature points 
on triangles see m) 

Adaptive Gaussian Quadrature. 

Special care is required close to EVs [8j, because the basis functions are no longer polynomials and 
the second order derivatives are singular at the EV. Eurthermore, the natural parametrization [28] 
produces only C'^-surfaces at EVs instead of C^-parametrizations. To obtain a -parametrization 
a suitable embedding has to be considered, e.g., using the characteristic map [Eli. For the imple¬ 
mentation, the evaluation algorithm of [ 2 H| can still be used, thanks to the integral transformation 
rule. 

Because of the structure of the subdivision scheme, the subdivision surface and correspondingly 
the basis functions <l>j are spline functions on local triangular mesh rings around each EV | 8 |. Thus, 
this subdivision ring structure can be used via an adaptive refinement strategy of the reference 
triangle A around the EV and an application of Gaussian quadrature on the resulting adaptive 
reference mesh to overcome the limitation of the standard Gaussian quadrature on the reference 
triangle. More explicitly, for triangles with an EV we decompose the associated reference triangle 
A into finer triangles Al {I = 1,... ,L) of level L (see Eigure left), perform the corresponding 
Gaussian quadrature GA(p) on all hner triangles A^ and call this adaptive Gaussian quadrature 
AG(p, L). Hence, the number of quadrature points K on these adaptively refined triangles depends 
on L: K = (3 • L -|- 1) • 12 for the Laplace-Beltrami equation 0 and K = (3 • L -|- 1) • 6 for the surface 
bi-Laplacian equation Q. It is worth mentioning that this expensive scheme has to be applied only 
for a small number of triangles around the finitely many EVs. Finally, let us remark that standard 
Romberg extrapolation does not lead to any improvement of the adaptive Gaussian quadrature 
due to the lack of smoothness of the integrands in the mesh size parameter and an adaption of the 
Romberg method for singular problems failed because the type of singularity is not explicitly known. 

Barycenter Quadrature. 

Implementation of the adaptive Gaussian quadrature is a tedious issue. In particular, in graphics 
where the achievable maximal order of consistency is not needed, already the computing cost of 
the usual Gaussian quadrature is significant. Therefore, reduced quadrature assembly is a common 
practice when implementing modeling or simulation tools based on NURBS as well as for subdivi¬ 
sion surfaces (e.g. mm)- The simplest and for subdivision surfaces widespread quadrature is the 
barycentric quadrature (BC) with the center point = ( 5 , 3 ) of the triangle A and the weight 
wi = ^. This rule is applied to regular as well as to irregular triangles and integrates only affine 
functions exactly. The method is also used in |T2| and leads to a reasonable consistency at least in 
the energy norm. 

Mid-edge Quadrature. 

An alternative, which shows superior performance with respect to the achievable convergence rates 
(cf. Section]^, is the mid-edge quadrature (ME) with K = 3 quadrature points at the midpoints 
of the edges, i.e., = (|, 0 ), = ( 5 , 2 ) and = ( 0 , g) with weights wi = W 2 = = g. 
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Figure 3: Left: decomposition of the unit triangle A in smaller triangles corresponding to the splines 
rings of the adaptive Gaussian quadrature AG(p, L) for L = 3. Right: illustration of the vertex 
configuration, with numbering as in Table (black straight and dashed lines) with a vertex pj of 
valence N in the neighborhood of the edge e connecting and P 2 - The red lines indicate the new 
vertices and edges in the 1-neighborhood of Pq after one level of subdivision. 


The ME quadrature integrates exactly polynomials of degree p = 2. Instead of following the direct 
implementation of (0, and ( |13[ ) we derive geometry-independent lookup tables of the 

basis function values and their derivatives on the unit triangle A at midpoints of the edges. This 
substantially simplifies the implementation. 

Figure (left) depicts a sketch of a generic local control mesh with two interior vertices, one 
possibly EV p^ of valence N and one regular vertex P 2 of valence 6. These two control points are 
surrounded by a fan of additional control points p^ {i = 3,N + A). To perform one step of 
subdivision, we introduce new vertices (which are always regular) on every edge and update the 
positions of the current vertices. The new vertex positions are linear combinations of the old vertices 
p^. Using the well-known position limit weights cj* for regular vertices [H] we can directly compute 
the limit position control vertex p^ in dependence of the coarse mesh vertices 

p^. This process applies to both the subdivision function values and their derivatives, which are 
listed in Table More explicitly, for the hat function Aj associated to a control vertex p^, the value 
of the basis function k) at the midpoint of the edge e £ le (connecting p^ and P 2 with local 

coordinates in the triangle (A, k) corresponding to the vertices Pi, P 2 , and pg) is a constant solely 
depending on N. The same holds for the derivatives ^i,i 2 and ^*,22 at in the 

directions on the reference triangle A. Let us remark that, if p^ and P 2 are EVs it is not clear if the 
global set of basis functions are linearly independent [25] . 

Based on Table the isogeometric subdivision approach can be implemented by iterating over 
all edges retrieving values from these lookup tables without implementing the box spline basis func¬ 
tions, their derivatives and the complex subdivision process itself. The iteration over edges instead 
of elements avoids to evaluate basis function values twice. Furthermore, the involved local matrices 
in the assembly of the mass and stiffness matrices are smaller for the mid-edge rule than for the 
barycenter rule. More explicitly, in the regular case, the local IgA-matrices have 144 entries for the 
BC rule versus 100 for the ME rule. Thus, based on the lookup tables, this leads to an overall faster 
assembly of mass and stiffness matrices than for the BC rule, even though the number of edges is 
I times the number of triangles on closed surfaces (e.g. see Table [^. Finally, let us remark that 
the mid-edge assembly process based on lookup tables can easily be generalized to other subdivision 
schemes like the Catmull-Clark 1201 or the Doo~Sabin scheme 
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<5^,1 (r,e) 

<^.,2(r,e) 


'^+I2(r,e) 

'^+22(r,e) 

i = 1 

eg-iQ-N-piN) 

192 

-19+16Af-/3CiV) 

24 

-19+16-iV-^C-^) 

48 

5-16-iV-/3(iV) 

4 

5-16-Af-^(iV) 

8 

-1 

i = 2 

62+16-/3(iV) 

192 

14-16-/3(Af) 

24 

14-16-/3(iV) 

48 

-2+16-/3(]V) 

4 

-l+8./9(iV) 

4 

-1 

i = 3 

25+16-/3CiV) 

192 

l-16-^(iV) 

24 

19-16-/3(iV) 

48 

-3+16-/3(]V) 

4 

-3+16-^CiV) 

8 

1 

2 

i = 4: 

2+16-^(Af) 

192 

-1-16-^CiV) 

24 

2-16-/3{Af) 

48 

4-p{N) 

-l+8./9(iV) 

4 

0 

i = 5, ■ • ■ , N—1 

16-/3(iV) 

192 

-2./3(iV) 

3 

-/3(JV) 

3 

4-p{N) 

2 • p{N) 

0 

i = N 

2+16-^(Af) 

192 

-l-16-^(iV) 

24 

-4-16-/3(Af) 

48 

4-p{N) 

l+8./9(iV) 

4 

1 

2 

i = AT + 1 

25+16-^(iV) 

192 

l-16-/3(iV) 

24 

-17-16-/3(Ai) 

48 

(-3+16./3(JV) 

4 

-3+16./3CJV) 

8 

1 

2 

i = N + 2 

3 

192 

1 

12 

-1 

48 

4 

-1 

8 

0 

i = N + 3 

1 

192 

1 

24 

1 

48 

4 

1 

8 

0 

i = N + 4 

3 

192 

1 

12 

5 

48 

4 

3 

8 

1 

2 

modification 
for = 3 



<^«,2(r,e) 


<^+i2(r,e) 

'^i,22iC, e) 

i = 3 

27+16-/3CiV) 

192 

-2/3(iV) 

3 

15-16-/3(iV) 

48 

-3+16./3(]V) 

4 

-l + 16.^(iV) 

8 

1 

i = 4: 

27+16-/3CiV) 

192 

-2/3(iV) 

3 

-15-16-/3(iV) 

48 

-3+16./3(iV) 

4 

-5+16./3(iV) 

8 

1 

2 


modification 
for A^ = 4 

MC,e) 


'^+2(r,e) 


'^+I2(r,e) 

^^i.22(C'’, e) 

i = 3 

25+16-/3CiV) 

192 

l-16-^(iV) 

24 

19-16-/3(iV) 

48 

-3+16./3(iV) 

4 

-3+16-^CiV) 

8 

1 

2 

i = 4 

4+16-^(Af) 

192 

-2-16-^(iV) 

24 

-2-16-^(N) 

48 

4-p{N) 

2 • p{N) 

1 

2 

i = 5 

25+16-^C^) 

192 

l-16-/3(iV) 

24 

-17-16-/3(Ai) 

48 

-3+16./3(]V) 

4 

-3+16./3(JV) 

8 

1 

2 


Table 1: Lookup tables for mid-edge quadrature rule ME. N denotes the valence of vertex and 
/3{N) = — (I + I cos For N = 3 and = 4 and some indices i the values differ 

from the those for general N. For the corresponding reference configuration see Figure (here i 
corresponds to p^). 


Remark Strang’s Lemma (see ga ia) provides the following error estimate for the numerical 
solution Uh of (14) 


\u - Uh\\ < 


\u - Uh\ 


+ 


\Uh - Uh\ 


Approximation error (Cea’s Lemma) Consistency error (Strang’s Lemma) 


< inf 
Vh&VS 


\u-Vh\\+ sup 


\a{vh,Wh) - a{vh,Wh)\ 


+ sup 


\e{wh) - iiwh)\ 




where u is the continuous solution of Q and Uh is the discrete solution of ([^. The last two terms 
measure the consistency of a and £f. Here, ||.|| := ||.||//i for the Laplace-Beltrami problem 0 and 
II-11 II-11/^2 for the surface bi-Laplacian problem Q. In fact, if a scheme with exact integration 
fulfills inf.u^gVh 11^ “ ''^hW < Ch^ (with p > 1 and a constant C), we ask for a numerical quadrature 
that preserves this order, i.e.. 


sup 

Wh&Vh 


\a{uh,Wh) - a{uh,Wh)\ \£{wh) - i{wh)\ j 


\Wh\ 




< ChP. 
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The optimal Wh = Wi’fi £ fixed h is given by Wh 

solution of 


, where Zh G VP is the 

y/a(Zh,Zh) 


a{zh, ^j) = a{uh, ipj) - a{uh, <Pj), V j G 

with Zhda = 0 and Zh = Figure 0 plots the resulting consistency error depending 

on the mesh size h of the control mesh. We observe an improved consistency by two orders for the 
mid-edge rule compared to the barycenter rule for the surface bi-Laplacian problem Q. 




log{ h) 


Figure 4: The consistency error is shown in a log-log plot for the Laplacian problem on the torus 
(cf. Fig. for varying grid size of the control mesh and for the barycenter quadrature (BC) and the 
midedge quadrature (ME). 

5 Numerical Applications 

We have implemented the proposed methods in C-|—|- and performed tests for four different subdi¬ 
vision surfaces: a torus surface with regular control mesh, a spherical surface (Spherical-3-4) with 
a control mesh with EVs of valence 3 and 4, a spherical surface (Spherical-5-12) with a control 
mesh with EVs of valence 5 and 12 and a complex real world hand model where the control mesh 
has altogether 119 EVs of valence 4, 5, 7, 8, 9 and 10. Let us emphasize that, in the spirit of the 
isogeometric approach, the control mesh determining the limit subdivision geometry is kept fixed. 
Eurthermore, the limit surface differs from a torus (Eigure with circular centerline and cross 
section or a perfect sphere (Eigure [^. Then, we successively refine the control mesh using subdivi¬ 
sion refinement to improve the accuracy of the discrete PDE solution. To run the simulations for 
Spherical-3-4, Spherical-5-12 and the hand model the initial control mesh has to be subdivided at 
least once to avoid EVs in direct neighborhood. In all tables, the mesh size h refers to the mesh size 
of these control meshes. 

On all of these surfaces we solve the two model problems Q and © for a given right-hand side / 
and compare the asymptotic error for different norms: the L^morm, the L7^-semi-norm and the 
semi-norm. To evaluate the experimental order of convergence (eoc), we have computed a reference 
solution for a control mesh with one additional level of global refinement compared to the finest 
mesh. Eurthermore to compute the reference solution we used the adaptive Gaussian quadrature 
AG{p,L) with L = 3 for Q and for L = 6 for Q (the estimated convergence rates did not change 
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Figure 5: Results for the torus: in the first row from left to right we depict the subdivision limit surface 
with control lines, isolines and color coding of the right hand side f{y) = sin(7ryi) sin(7r7/2) sin(7ry3) 
(—1.0 ■ 1.0), the (numerical reference) solution of the Laplace-Beltrami problem (—0.054 1 ■ 

0.054), and the (numerical reference) solution of the surface bi-Laplacian problem (—0.003 ; ■ 

0.003) on A4; in the second and third row log-log plots of the error are reported for the Laplace- 
Beltrami problem (2’*'^ row) and the surface bi-Laplacian problem (3'"'^ row), respectively (from left 
to right: L^-norm, Lf^-semi-norm and Lf^-semi-norm). 


for larger L). As a consequence, the error plots for the finest discretization are less reliable as it 
becomes apparent in Figure]^ for the L^-error plot for problem Q and the spherical shape with low 
valence EVs (Spherical-3-4). For all the experiments, double precision arithmetic is used. 

Figure shows our results for the torus with regular control mesh without EVs. All quadrature 
rules achieve optimal convergence rates (cf. [50]) for the second-order problem Q, i.e., 4 in the L^- 
norm, 3 in the iL^-semi-norm and 2 in the Lf^-semi-norm (with non-adaptive Gaussian quadrature 
GA(p)). Here, ’’optimal” reflects what we expect for the quartic box spline |5D|. Eor the fourth-order 
problem Q all quadrature rules achieve optimal rates 4 in the L^-norm, 3 in the Lf^-semi-norm and 
2 in the H^-semi-norm (with non-adaptive Gaussian quadrature GA(p)), except the BG rule. Let us 
remark that the C^-regularity allows for an Aubin-Nitzsche type argument |44j to obtain estimates 
in norms other than the energy norm. In all cases GA(p) performed best. 

The numerical results of the spherical control mesh with EVs of valence 3 and 4 (Spherical-3- 
4) depicted in Figure]^ For the second-order PDE Q Gaussian quadrature and adaptive Gaussian 
quadrature achieve the optimal order of convergence. On finer resolutions ME and BG do not achieve 
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log( h ) log( h ) 


Figure 6; Results for Spherical-3-4; As in Figure we show in the top row the limit surface, 
f{y) = sin(37ryi) sin(37r?/2) sin(37ry3) ( — 1.0 ; O 1.0), the solution of the Laplace-Beltrami problem 
(—5.25e —3 ; a 5.25e —3), and the solution of the surface bi-Laplacian problem (—9.78e —5 I 3 
9.78e — 5), together with the error-plots for the Laplace-Beltrami problem (2"’'^ row) and the surface 
bi-Laplacian problem (3^*^ row) (from left to right: L^-norm, 77^-semi-norm and 77^-semi-norm). 


optimal convergence rates. GA, AG and ME achieve similar behavior for the fourth-order problem 
and again the BC method performs worse. The reduced convergence rate in the L^-norm for the 
fourth-order problem is expected to reflect the reduced regularity H 77^, which precludes the 

application of an Aubin-Nitzsche type argument. 

For the control mesh with EVs of valence 5 and 12 (Spherical-5-12) the convergence results are 
depicted in Eigurej^and coincide with our general findings for control meshes with EVs of valence 
greater than 6. Here, all quadrature schemes show a similar performance, for the Laplace-Beltrami 
problem Q: 3 in the L^-norm, 2 in the 77^-semi-norm and 1 in the 77^-semi-norm, and for the 
surface bi-Laplacian problem ([^ : 2 in L^-norm, 2 in the 77^-semi-norm and 1 in the 77^-semi-norm. 
This observed loss of approximation order compared to the quartic box spline coincides with the 
theoretical work by m with the reasoning that Loop subdivision functions cannot reproduce cubic 
polynomials around extraordinary vertices of valence greater then 6. 

Finally, we consider the hand model as a complex subdivision surface with many EVs in Figure 
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Figure 7: Results for Spherical-5-12: We follow the presentation in Fig. [^and plot the limit surface 
with IgA-element lines, the (same) right hand side function / now on this surface, the solution of the 
Laplace-Beltrami problem (—1.58e — 2 ; 31 1.58e — 2), and the solution of the surface bi-Laplacian 

problem (—1.27e — 3 ; 31 1.27e — 3), again together with the error-plots for the Laplace-Beltrami 

problem (2"''^ row) and the surface bi-Laplacian problem (3^“^ row) (from left to right: L^-norm, 
77^-semi-norm and 77^-semi-norm). 

Given the complexity of the geometry and the number and valence of the extraordinary vertices both 
for the Laplace-Beltrami and surface bi-Laplacian problem the asymptotic regime of error reduction 
seems not to be reached even though the reference solution is computed on a mesh with 750A: vertices. 
The BC fails for Bi-Laplacian probably caused by the increased number of EVs and their partially 
high valences. Figure shows eigenfunctions of the Laplace-Beltrami operator computed via inverse 
vector iteration with projection. The depicted results underline that the methods discussed so far 
are also applicable in the numerical eigenmode analysis, which turned out to be an indispensable 
tool in geometric data analysis and modeling. 

Compared to Gaussian quadrature and in particular to the adaptive Gaussian quadrature the 
mid-edge and the barycenter quadrature are significantly cheaper as reported in Table Because 
of our implementation, the mid-edge rule based on lookup tables and assembly via an edge-iterator 
performs even better than the (non-optimized) barycenter rule. 


14 




































































Geometry 

\^v 1 

GA(p) 

AG{p,L) 

— ^M { — ^mY 

— ^M 

BG 

— ^M 

ME 

Spherical-5-12 

9218 

1.209 

1.417 

1.374 

1.563 

0.164 

0.272 

0.107 

0.136 

Hand 

11586 

1.914 

2.233 

2.615 

3.424 

0.269 

0.460 

0.157 

0.198 


Table 2: Comparison of assembly times for the stiffness matrices (10) and (11) and the different 


quadrature rules (Gaussian GA{p), adaptive Gaussian AG{p,L), barycenter (BC) and mid-edge 


(ME)). All times reported in seconds. 









Figure 8: Results for complex hand shape: We follow the presentation in Fig. and plot the limit 
surface with IgA-element lines, the (same) right hand side function / now on this surface, the 
solution of the Laplace-Beltrami problem (—1.47e — 2 ' 3 1.14e — 2) and the solution of the 

surface bi-Laplacian problem (—8.14e.4 : 3 7.56e.4), again together with the error-plots for the 

Laplace-Beltrami problem (2”'^ row) and the surface bi-Laplacian problem (3'’’^ row) (from left to 
right: L^-norm, L7^-semi-norm and iL^-semi-norm). 
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6 Conclusions 

We have explored the discretizations of geometric partial differential equations on subdivision sur¬ 
faces by means of an isogeometric approach. To this end, we have focused on Loop’s subdivision 
method and have investigated the impact of different numerical quadrature rules on the expected 
convergence rates. As expected, there is a trade-off between robustness and computational effort. 
The adaptive Gaussian quadrature turns out to be the most robust in the carefully selected set of 
test cases but at the same time computationally demanding for instance for real-time applications in 
modeling and animation. In fact, the mid-edge quadrature based on lookup tables has been singled 
out as a promising compromise for both graphics and engineering applications. It is very easy to 
implement and it performed well in all considered test cases. Fourth-order problems like the surface 
bi-Laplacian problem considered here, are by far more critical with respect to the accuracy of the 
numerical approximation than second-order problems. In general, adaptivity around EVs is essential 
to ensure the overall robustness on very fine computational meshes. 
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